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Abstract 

The modified mixing length (MML) turbulence 
model was installed in the Proteus Navier-Stokes 
code, then modified to make it applicable to a wider 
range of flows typical of aerospace propulsion appli- 
cations. The modifications are based on experimental 
data for three flat-plate flows having zero, mild 
adverse, and strong adverse pressure gradients. 

Three transonic diffuser test cases were run with the 
new version of the model in order to evaluate its per- 
formance. All results are compared with experimen- 
tal data and show improvements over calculations 
made using the Baldwin-Lomax turbulence model, 
the standard algebraic model in Proteus. 

Nomenclature 

A + van Driest damping constant = 26 
Cf local skin friction coefficient 

C x MML parameter; controls mixing length 

saturation level 

C 2 MML parameter; controls curvature of 

blending region 

Baldwin-Lomax turbulence model constant 

cp 

= 1.6 

c Kleb Baldwin-Lomax turbulence model constant 

= 0.3 

C wk Baldwin-Lomax turbulence model constant 

= 0- 25 

Dj, Dj parameters used in turbulence model 

averaging for multiple walls 
fj, f 2 parameters used in turbulence model 

averaging for multiple walls 
F(y) function in Baldwin-Lomax turbulence 

model 

F Kleb Klebanoff intermittency factor 

F max parameter in Baldwin-Lomax turbulence 

model 
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parameter in Baldwin-Lomax turbulence 
model 

MMLPG parameters 

throat height of Sajben diffuser 

coefficient of thermal conductivity 

turbulent mixing length 

static pressure 

total pressure 

vector of dependent variables 

ratio of exit static pressure to inlet total 

pressure for Sajben diffuser 

Reynolds number based on x-coordinate 

time 

velocities 

fireestream x-velocity 
shear velocity 
total velocity 

difference between maximum and 

minimum total velocities 

Cartesian coordinates 

y coordinate nondimens ionalized by shear 

length scale 

shear length scale 

parameter in Baldwin-Lomax turbulence 


model 

P Clauser’s equilibrium parameter 

5 boundary layer thickness 

S] displacement thickness 

t£^\ ££^second- and fourth-order artificial 

viscosity coefficients in constant coeffi- 


cient model 

Zj implicit artificial viscosity coefficient 

k von Karman constant = 0.4 

k 2 , K 4 constants in nonlinear coefficient artificial 
viscosity model 

X second coefficient of viscosity 

\i molecular viscosity 

tj computational coordinate directions 
p density 


1 


c 

pressure gradient scaling parameter in 
nonlinear coefficient artificial viscosity model 

V 

spectral radius in nonlinear coefficient artifi- 
cial viscosity model 

X 

shear stress 


shear stress at interior grid points 

CO 

vorticity 

Subscripts 

cap 

capping or saturation value 

e 

edge of boundary layer 

eff 

effective 

i,j 

indexes in the x and y directions 

inner 

inner region of boundary layer 

max 

maximum 

min 

minimum 

outer 

outer region of boundary layer 

t 

turbulent 

w 

wall 

Superscripts 

+ 

nondimensionalized by the shear length scale 
Introduction 


The Proteus Navier-Stokes code 1 was developed 
for aerospace propulsion flow applications and solves 
the Reynolds-averaged, unsteady, compressible Navier- 
Stokes equations in strong conservation-law form. It 
uses a fully-coupled alternating direction implicit 
solution procedure 2 with generalized first or second 
order time differencing, 3 employs implicit boundary 
conditions, and linearizes all terms using second order 
Taylor series expansions. Two turbulence models are 
available in Proteus: The Baldwin-Lomax algebraic 
model (BLM) 4 and the Chien k-e model. 5 

For engineering applications, algebraic models 
offer the most algorithmically simple and computation- 
ally efficient approach to turbulence modeling. Though 
two-equation models are sometimes thought to be more 
accurate at calculating complex flows, they are more 
difficult to implement, require initial and boundary 
conditions that are often difficult to define, and are 
usually more computationally expensive than algebraic 
models. The Baldwin-Lomax model (presented in 
Appendix A) is the most widely used algebraic model, 
even though it is known to have difficulties computing 
flows with strong adverse pressure gradients and large 
separated regions. 6 ” 14 The modified mixing length 
model (MML) was developed specifically to handle the 
separation that occurs on airfoils with leading edge ice 
accretions, and it produced significantly better results 
than BLM for such flows. 6 The success of these calcula- 


tions warrants further evaluation and development of 
MML. 

The objective of this work is to evaluate and then 
modify MML to improve its performance for adverse 
pressure gradient flows. To accomplish this, MML was 
installed in Proteus and first used to calculate zero 
pressure gradient flow over a flat plate. It was modified 
to more accurately predict the boundary layer growth for 
this flow; then the experimental data of Bradshaw 15 was 
used to make additional modifications for adverse 
pressure gradient flows. The resulting model, called 
MMLPG, was then validated for three transonic difiFuser 
test cases. The CRAY Y-MP computer at the NASA 
Lewis Research Center was used for all calculations. 

The MML Model 

The modified mixing length model was devel- 
oped by Potapczuk 6 to fill the need for an algebraic 
model to calculate turbulent flow with large separated 
regions. The particular problem of interest was an airfoil 
at angle of attack with and without leading edge ice 
accretions. Previous calculations made with BLM 
produced nonphysical results; the source of the problem 
was the function F(y) (defined in Appendix A), which 
could have multiple peaks for this flow. This behavior is 
shown in figure 1, taken from reference 6. As the 
relative magnitudes of the F(y) local maxima change, 
y max may suddenly jump, producing unrealistic disconti- 
nuities in the turbulent viscosity. Selection of the global 
maxima often results in a gross over-prediction of the 
turbulent viscosity. Some authors®’ ® have found that 
choosing the outermost peak produces better results, 
while others have elected to use the innermost peak. 10 

The MML model avoids the need to seek a 
maximum of some ad hoc function. In accordance with 
Prandtl’s mixing length theory, the MML model deter- 
mines the mixing length using the wall shear stress and 
the normal distance from the wall, with the maximum 
mixing length capped at a given value. As a result, it is a 
two layer model; the length scale depends on conditions 
near the surface and remains constant in the separated 
region. This assumption is valid since there is no 
substantial enhancement of turbulence in separated 
regions. The turbulent viscosity, is given by 

n, - piVi w 

where p is the density, / is the mixing length and |co| is 
the vorticity magnitude. Figure 2, taken from reference 
16, shows the behavior of the mixing length in a turbu- 
lent boundary layer. Several empirical formulas are 
available to evaluate the inner region, 16, 17 which 
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consists of the viscous sublayer and the overlap layer 
The MML model uses the van Driest formulation, 
given by 


= *y 
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( 2 ) 


where A + = 26 and the value of k, the von Karman 
constant, is 0.4. The quantity y + is defined as 



where y* is the shear length scale 
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very large. To avoid this problem, the following local 
average was incorporated: 6 

N = 01 l x i-2l +0 - 2 IV.I +0 - 4 N +0 - 2 l T i + tl +0 - 1 l t i + 2|( 8 ) 

The subscripts in equation (8) refer to grid points in the 
streamwise direction along the wall. 

Calculations for iced airfoil flows made by 
Potapczuk with the MML model showed improvements 
over BLM calculations for predictions of the separated 
region, the maximum lift coefficient and vortex shedding 
frequencies. Since the MML model was developed to 
solve the specific problem of flow over airfoils, a 
comprehensive evaluation of the model for more general 
flow fields was not a part of that study. The objective of 
the present study is to evaluate the MML model for 
general zero pressure gradient and adverse pressure 
gradient turbulent boundary layer flows and examine 
possible modifications to improve the performance of the 
model. 


Here, p is the molecular viscosity and x w is the wall 
shear stress. For y + £ 5 A + (but still in the “inner” 
region), the mixing length is approximated by Ky ; this is 
the original Prandtl theory, and is consistent with the 
well-known logarithmic profile. In the outer region of 
the boundary layer, the outer mixing length is assumed to 
behave according to 


/outer = constant Xy* 


(5) 


The MML model uses a blending function to 
give a smooth transition between the inner and outer 
layers; this is given by 


/(y> = 



r y +> 

( f + 0^ 

+ 

1 1 _ 1 1 _ y_ i 

i- c A 

l l 1 C,J J 


^ > 

It 

.ftP 

T* 

y + 2C, 


y + <C| (6) 


(7) 


In this formulation, C^y* is the distance from the surface 
at which / saturates, and C 2 controls the curvature of the 
blending region. Figure 3 shows a typical MML model 
mixing length profile for an attached wall boundary 
layer. 

Near separated regions, x w approaches zero 
which causes y* and thus the mixing length to become 


Evaluation and Modi fication of MML 
Evaluation of MML for Zero Pressure Gradient Floai 

The test case of incompressible, zero pressure 
gradient, turbulent flow over a flat plate, as shown in 
figure 4, was used to evaluate MML. The grid, shown in 
figure 5, had 51 points in both the streamwise and 
normal directions and had grid points clustered at the 
wall to resolve the boundary layer and at the upstream 
boundary to resolve the imposed boundary condition. In 
addition, it was evaluated to insure grid independence 
for zero pressure gradient flow. The reference velocity, 
temperature, pressure and length used in Proteus were 
33.53 m/s, 288.3 K, 101.3 kPa and 1.98 m, respectively. 
At the upstream boundary, the velocity profile, which 
was computed using the correlation of Musker,^ was 
held fixed The flow was computed using both MML 
and BLM. The MML constants were chosen as (^=3000 
and C 2 = 5 , which were found to give good results at 
Re x =7xl0 6 . Both turbulence models produced good 
agreement with experimental velocity-defect profiles,^ 
as shown in figure 6 at Re x =7xl0 6 . The quan tity u T in 
figure 6 is the shear velocity, given by u t *= ^(jtj/p). 
At other locations on the plate (i.e., at other Reynolds 
numbers) the BLM velocity-defect profiles correctly 
exhibit similarity but the MML profiles do not, as shown 
in figure 7. 

In a turbulent flow over a flat plate, the boundary 
layer thickness increases with increasing x-distance 
along the plate. To accurately model this flow, the turbu- 
lent length scale must also increase proportionately with 
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the boundary layer thickness. In MML t the outer length 
scale, as given in equation (7) is equal to a constant 
times the shear length scale, y* . The increase in y* with 
x-distance is almost negligible, resulting in an essen- 
tially constant value of the outer length scale for all 
Reynolds numbers. Plots of p t , as shown in figure 8, 
illustrate that the turbulent viscosity profiles calculated 
with MML are nearly the same at all Reynolds numbers, 
but the BLM profiles increase with increasing 
Reynolds number. Though MML produced the correct 
length scales for an airfoil near stall, 6 modifications are 
needed to make it applicable to general boundary layer 
flows. 

In order to make MML applicable over a range of 
Reynolds numbers, the optimal saturation lengths, or C\ 
values, were found at several Reynolds numbers. The 
following simplified formulas were used to calculate the 
inner and outer mixing lengths: 


l + = Ky + 


f y + ^ 


1 -c 


c. <y 


I cip — kCj, C^y 


(9) 

( 10 ) 


Here, / + is the nondimensional form of the mixing 
length, equivalent to l/y* , and the outer length scale, 
/ + cip , is simply the inner length scale evaluated at 
y + = C, . From these results, / + cip was found as a 
function of the skin friction, Cf, giving 


/ + cp = 1860- 6.20 xl0 5 c f (11) 

The velocity-defect profiles of figure 9 show that 
equations (9) and (11) with / + = rain (/*, / + C « P ) . allow 
the mixing length to grow proportionately with the 
boundary layer thickness. The modified MML is better 
than BLM at predicting the local skin friction coefficient, 
Cf, as shown in figure 10; the wiggles at the upstream 
boundary are a result of the imposed upstream boundary 
condition. 

Modifications for Adverse Pressure Gradient Flows 
Two equilibrium pressure gradient flows of 
Bradshaw 15 were used to modify MML for adverse 
pressure gradients effects. Equilibrium turbulent flows 
are flows which have a constant value of Clauser’s 
equilibrium parameter, 21 


B - !i- e 


( 12 ) 


In addition, they correspond to a power-law velocity 
profile distribution, u c « x\ where the magnitude of the 
exponent, a, indicates the strength of the pressure gradi- 
ent. They also exhibit similarity when plotted in veloc- 
ity-defect coordinates. Three flows were examined in 
the experimental study of Bradshaw; 15 these were flows 
with zero, mild, and strong adverse pressure gradients. 
The corresponding values of the exponent, a, are 0, 
-0.15, and -0.255, respectively; the corresponding values 
of P are 0, 1 and 5. 

The modifications to the turbulence model are 
based on the trends exhibited in the mixing length at the 
three pressure gradients as shown in figure 1 1 taken from 
Bradshaw. 15 Note that for all three pressure gradients, 
the maximum mixing length is approximately 0.088, the 
saturation distance from the wall is approximately 0.45, 
and the slope of the curves near the wall increases with 
the strength of the pressure gradient. These three 
features were used to develop the following model: 



y + <G, (13) 


^G 1, y >G > 


(14) 


A new parameter, G3, has been introduced and the 
constants Cj and Cj in the original MML model have 
been replaced by the functions Gj and 62, where 


G, = 0.4G 4 

(15) 

G 2 = 5G 3 k 

(16) 


Here, G 4 is essentially a nondimensional boundary layer 
thickness which is a function of P and Cf, and G3 controls 
the slope of the mixing length curve and is a function of 
p. The following correlation was assumed for G 4 : 


G 4 - G 5 + G 6 c f 


(17) 


Separate values of the parameters G3, G5, and Gg, corre- 
sponding to each of the three pressure gradients, were 
found and are given in table 1. This results in essentially 
three separate models, one for each pressure gradient, 
depending on the set of parameters used. The results of 
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these modifications are compared with Baldwin-Lomax 
calculations in the velocity-defect plots of figures 12 
through 14. (Note: The Baldwin-Lomax results for the 
zero pressure gradient case are given in figure 7.) The 
reference conditions used are the same as those given 
earlier for the zero pressure gradient case. For the two 
adverse pressure gradient cases, the turbulent velocity 
profiles at the upstream boundary were computed using a 
cubic spline fit of the Bradshaw experimental data and 
held fixed; the appropriate pressure gradient was 
imposed at the freestream boundary. Both cases were 
computed using the grid of figure 5, but for the strong 
pressure gradient case, the number of grid points in the 
vertical direction was increased to 101. 

Final Model 

The final step in developing this turbulence 
model was to combine all of the above modifications into 
one general turbulence model. To accomplish this, the 
following correlations were developed for the parame- 
ters G3, G 5 and Gg: 


g 3 = 1.0 

0<O.O 

(18.a) 

G3 = 1.0 + 0.307 P- 0.039 lp 

2 0.05 055.34 

(18.b) 

G 3 = 1.52 

0 > 5.34 

(18.c) 

G 5 = 23, 300 

0<O.O 

(18.d) 

G s = 23, 300 + 85600 - 123O0 2 0.0 5 0 5 5.34 

(18.e) 

G 3 = 33,900 

0 > 5.34 

(18.0 

G 6 = -7.75xl0 6 

0<O.O 

(18.g) 

G 6 = - 7.75 x 10 6 - 4.5 1 x 1O 6 0 + 386, OOO0 2 

0.0 £05 5.34 

(18 .h) 

G 6 = -20, 900 

0 > 5.34 

(18.i) 


The available experimental data are limited to only the 
three values of P which are in the range 0 £ P £ 5.34, and 
the quadratic correlations of equations (18.b), (18.e) and 
(18.h) are based on this limited data. For P<0, the values 
in (18.a), (18.d) and (18.g), were obtained by evaluating 
the quadratic equations at P=0. Similarly, for P>5.34, 
the values in (18x), (18.0 and (18.i) were obtained by 
evaluating the quadratic equations at P=5.34. 

Since P, defined in equation (12), is a function of 
the displacement thickness, 5j, a correlation was also 


developed to avoid the problem of calculating 5| directly 
and thus having to define the edge of the boundary layer. 
This resulted in 

6, = (G ? + G g c f ) y* (19) 

The parameters G7 and Gg were defined in a manner 
similar to G3, G5 and G^ as given below. 


G ? = 2910 0<O 

(20.a) 

G 7 = 2910 + 27000 - 3430 2 0 5 0 5 5.34 

(20.b) 

G ? = 7560 0 2: 5.34 

(20.c) 

G g = -96900 0 < 0 

(20.d) 

G g = - 988, 000 - 1.15 x 1O 6 0 + 89, OOO0 2 

05055.34 

(20.c) 

G g = —4.57 x 10 6 0 > 5.34 

(20.0 


The value of P used to define G7 and Gg is lagged in 
time. 

The final model with pressure gradient modifica- 
tions, called MMLPG, was developed using the equilib- 
rium turbulent flows of Bradshaw and is defined by 
equations (13) through (20). The resulting velocity- 
defect profiles for all three pressure gradient flows are 
shown in figure 15 and exhibit good agreement with the 
experimental data, with the exception of the strong 
pressure gradient case. The calculations were performed 
on a CRAY Y-MP computer and the computational times 
are given in table 2. The strong pressure gradient case 
took considerably longer to reach convergence because 
the code had difficulties resolving oscillations induced at 
the upstream boundary, which used a fixed velocity 
profile for the inflow boundary condition. 

Averaging for Multiple Boundaries 

If both walls in a given coordinate direction are 
solid surfaces, the turbulent mixing lengths are 
computed separately at each surface and then averaged. 
The Sajben diffuser, which is described in the next 
section, has solid walls at the upper and lower vertical 
boundaries, and is a typical example of a geometry 
which would require averaging of the mixing length. 

The averaging formula of Appendix A, equation (A. 10), 
which was used to average the F waicc function in the 
Baldwin-Lomax model, is also used here to average the 
mixing length: 
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/ = 


l l ( [ +l 2 f 2 

f i+ f 2 


( 21 ) 


If the lower and upper boundaries in the vertical direc- 
tion, are solid surfaces, as in the Sajben diffuser, then l\ 
and / 2 are the mixing lengths at the lower and upper 
boundaries, respectively. The functions fj and f 2 are 
defined in equation (A. 11) of Appendix A. 


Adverse Pressure Gradient Test Case s 
To evaluate MMLPG for some typical propulsion 
flows, a converging-diverging duct, referred to as the 
Sajben diffuser, was used. This duct was designed to be 
representative of the diffuser portion of the inlet for a 
rocket/ramjet propulsion system. Detailed experimental 
and computational data are available for flows with and 
without external excitations. 22 " 27 This study, however, 
dealt only with the unexcited flows. The geometry of the 
diffuser is given in figure 16: the throat height, H, is 44 
mm; the entrance-to-throat height ratio is 1.4, and the 
exit-to-throat height ratio is 1.5. The grid, shown in 
figure 17, is the same as that used by references 1 and 26, 
and has 81 streamwise points and 51 vertical points. It 
was packed in the vertical direction near the walls in 
order to resolve the turbulent boundary layers and in the 
streamwise direction near the throat to resolve the strong 
gradients. The reliability of this grid is discussed in 
Appendix B. Three transonic flow cases were run. The 
flowfields were determined by setting R, the ratio of the 
exit static pressure to the inlet total pressure. The first 
case had a weak normal shock with R=0.82; the second 
case had subsonic flow throughout (no shock) with 
r= 0.862, and the third case had a strong normal shock 
with R=0.72. The reference velocity, temperature, 
pressure and length used in Proteus were 4.72 m/s, 292 
K, 135 kPa, and .044 m respectively. These values 
match the values used in other numerical simulations of 
this flow. 1, 24f 26 The initial conditions were zero velocity 
and constant temperature and pressure everywhere in the 
flow field. Both cases were run using MMLPG and the 
Baldwin-Lomax model. 


Weak Shock Case 

The weak shock case was used as an example 
case in the Proteus User’s Manual, 1 and therefore was 
run first in order to gain familiarity with running this 
type of flow. It was computed as described in reference 
1: First the exit pressure was gradually reduced to R = 
0.1338 to establish supersonic flow throughout the 
diffuser; then it was gradually raised to R = 0.82, the 
desired ratio to establish the weak normal shock, and 
iterated until the solution was no longer changing appre- 


ciably with time. A plot of the static pressure on the top 
wall at two locations, one upstream and one downstream 
of the normal shock, as the solution progresses is shown 
in figure 18. This indicates that pressure reaches a 
steady state level, which, for practical engineering 
purposes, can be considered a converged solution. A 
closer examination of the results indicates that the 
solution oscillates slightly about a mean steady level. 
This may be caused by inherent unsteadiness in the flow; 
Salmon et al. 22 mention that very low-amplitude, self- 
sustaining oscillations were observed experimentally. It 
is more likely, however, that the oscillations present in 
this calculation are numerical in nature, which is 
common for flows with shock waves. The oscillations 
originating at the shock may not be damped out by the 
artificial viscosity and therefore tend to migrate 
upstream. The artificial viscosity used in Proteus to 
calculate this flow was second- and fourth-order explicit, 
both using the nonlinear coefficient model of Jameson et 
al.; 2 ® the respective smoothing coefficients are k 2 and 
K4, as given in Appendix B. For the entire calculation, 
k 2 was set to 0.1; K4 was set to .005 for the first 6000 
iterations, while the exit pressure was changing, and 
decreased to .0004 for the remaining 3000 iterations, 
which were at a constant exit pressure. More details 
about the effects of the artificial viscosity on this 
solution are provided in Appendix B. 

The static pressure distribution on the top and 
bottom walls is given in figure 19. The shock location 
on the upper wall and the shock Mach number at the 
edge of the upper wall boundary layer are given in table 
3. Both MMLPG and BLM accurately predict the 
pressure distribution on the wall and the location of the 
shock. Each case was run for 9000 iterations and calcu- 
lations made using MMLPG and BLM required 
3.44xl0’ 5 sec/iteration/grid point and 3.36x1 0'^sec/itera- 
tion/grid point respectively on the CRAY Y-MP 
computer. 

No Shock Case 

The second diffuser test case did not have a 
normal shock wave. To compute this case, the exit 
pressure was gradually lowered to R =0.862 then iterated 
until the solution stopped changing. A steady state 
solution was reached with subsonic flow throughout the 
entire diffuser. The Proteus default artificial viscositv, 
which uses the constant coefficient model of Steger 29 
with both fourth-order explicit and second-order implicit 
artificial viscosity, was used; the smoothing coefficients, 
E£^ and £/ (defined in Appendix B), had values of 1.0 
and 2.0, respectively. 
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The static pressure distribution on the top and 
bottom walls is shown in figure 20 and indicates that 
MMLPG is clearly better than BLM at predicting the 
pressure distribution, though it still predicts a larger 
pressure drop than that indicated by the experimental 
data. The MMLPG results are similar to the calculations 
of Hsieh et al 25 who attributed the lower throat pressure 
to the fact that the experiment was highly sensitive to 
small perturbations in exit pressure. The maximum 
Mach numbers in the diffuser are given in table 4. 
Although no experimental data is available to compare 
these values, the MMLPG results are in best agreement 
with the calculations of Georgiadis 27 who used the 
PARC Navier-Stokes code 30 for the same geometry. Of 
the three diffuser test cases, this flow is most similar to 
the benchmark cases used to derive MMLPG, although 
its pressure gradients are stronger, and indicates that 
MMLPG is capable of computing the flows for which it 
was designed. Each case was run for 9000 iterations and 
calculations using MMLPG and BLM required 3.45x10 3 
sec/iteration/grid point and 3.36x10"'* sec/iteration/grid 
point, respectively, on the CRAY Y-MP computer. 

Strong Shock Case 

The final diffuser flow computed was the case 
with a strong normal shock positioned in the throat. This 
case was run in a manner similar to the weak shock case: 
First the exit pressure was gradually lowered to 
R=0.1338 to achieve supersonic flow throughout the 
diffuser; then it was gradually raised to R=0.72 to estab- 
lish the strong normal shock in the throat region, and 
iterated there until the solution stopped changing appre- 
ciably with time. Proteus was run in both steady and 
unsteady modes to try to simulate the experimentally 
observed self-excited oscillations of 217 Hz. 22 Unsteady 
mode in Proteus is achieved by calculating a global time 
step whereas steady mode uses a local time step to speed 
up the computation. Neither calculation simulated the 
experimentally observed oscillatory behavior, but 
instead produced very small numerical oscillations in the 
flow properties. (The artificial viscosity used for the 
strong shock calculations was the same as that used for 
the weak shock calculation.) Figure 21 shows the static 
pressure on the top wall at the experimental shock 
location and illustrates the behavior of these small oscil- 
lations; the calculation shown was run in unsteady mode 
with MMLPG. 

The static pressure on the top and bottom walls 
are presented in figure 22 and the shock location and 
Mach number at the edge of the top wall boundary layer 
are given in table 5. Both MMLPG and BLM predicted 
the shock location too far downstream. The experiment 


predicted a region of separation on the top wall just 
downstream of the shock with the flow reattaching at xJ 
H=6.0. MMLPG predicted a very small region of 
separation on the top wall which reattached at x/H=3.6. 
BLM predicted very small regions of separation on both 
the top and bottom walls which reattached at x/H=3.8 
and x/H=6.2. The separated behavior is illustrated in 
figure 23, which gives the velocity profiles at four 
locations downstream of the shock. These peculiar 
results can be attributed to the fact that both MMLPG 
and BLM compute very high values of m due to the large 
increase in vorticity downstream of the shock. The poor 
performance of both models for this case can also be 
attributed to the fact that all of the models are equilib- 
rium turbulence being used to calculate a flow which is 
clearly nonequilibrium. The inadequate performance of 
MMLPG for the strong shock case is also explained by 
the derivation of the model, which is based on experi- 
mental data for beta values between 0 and 5, while this 
flow encountered P values as high as 12,000. Each case 
was run for 10,000 iterations and the steady calculations 
using MMLPG and BLM required 3.54xl0" 5 sec/itera- 
tion/grid point and 3.82x10’* sec/iteration/grid point, 
respectively, on the CRAY Y-MP computer. 

Summary and Conclusions 
In the current work, the range of applicability 
of the MML algebraic turbulence model was extended to 
calculate more general boundary layer flows with zero 
and adverse pressure gradients. To accomplish this 
objective, MML was first modified to predict the appro- 
priate boundary layer growth with increasing Re x by 
incorporating a relationship for the mixing length as a 
function of the local skin friction coefficient. Using the 
experimental data of Bradshaw for zero, mild and strong 
adverse pressure gradient flows, modifications were also 
added to account for adverse pressure gradient effects. 
The resulting generalized model, called MMLPG, 
accurately predicted zero and adverse pressure gradient 
flows over a plate and exhibited better agreement with 
experimental data than the Baldwin-Lomax model. 

To more thoroughly evaluate MMLPG for other 
adverse pressure gradient flows, this model was also 
used to calculate three transonic dtffuser flow test cases: 
flow with a weak shock, flow with no shock, and flow 
with a strong shock. For the weak-shock case, MMLPG 
and BLM did equally well in predicting the shock Mach 
number and location, and also in predicting the static 
pressure distribution on the top and bottom diffuser 
walls. For the no shock case, MMLPG was significantly 
better than the Baldwin-Lomax model at predicting the 
static pressures on the walls and at predicting the 


7 



maximum Mach number in the duct. Both models inade- 
quately predicted the strong shock flow, over-predicting 
the shock Mach number and location and under predict- 
ing the size of the separation on the top wall. In 
addition, BLM predicted a small separation on the 
bottom wall, although no separation was observed there 
experimentally. This strong shock flow over-stepped the 
bounds of the assumptions made in the development of 
both of these equilibrium turbulence models. 

Overall, the flat plate and transonic diffuser 
results indicate that MMLPG is capable of accurately 
predicting turbulent flows with and without adverse 
pressure gradients. Future work should include contin- 
ued validation of the model for these types of flows as 
well as continued development of the model to better 
account for stronger adverse pressure gradient flows both 
with and without separation. 
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A ppendices 

A ppendix A: The Baldwin-Lomax Turbulence Model 
A generalized version of the Baldwin-Lomax 
algebraic turbulence model 4 is available in Proteus. 1 
The turbulent shear and normal stresses and the turbulent 
heat flux are modeled using the Boussinesq approach, 
where the effective viscosity is defined as p eff = p + p t , 
the second coefficient of viscosity is defined as 
X tff = X + X t , and the effective thermal conductivity 
coefficient is defined as k cff = k + k t . 

For wall bounded flows, the Baldwin-Lomax 
model is a two-layer model: 


* 

y ^ ycrouover 

Pouter ■ 

y > ycrosfover 


where y^^over * s smallest value of y at which the inner 
and outer region values of p t are equal. For free turbu- 

lent flows, n t = (Ji,) 0Bter . 

1 . Inner Region 

The inner region turbulent viscosity is computed 

from 


<H,>. 


p/ 2 |to) 


where / is the mixing length given by 


+ \ 


/ = Kyi 


1-e 


(A2) 


(A3) 


The quantity |o)[ is the magnitude of the total vorticity, 
defined for two-dimensional planar flow as 
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N 


dv 9uj 
9x 9yl 


(A.4) 


2. Outer Region 

In the outer region, the turbulent viscosity is 
given by 


^ outer ^ C cpP F Klcb F wakc 

where K is the Clauser constant, taken as 0.0168 and C cp 
is a constant taken as 1.6. The quantity F wakc is 
computed from 


vertical direction, if the upper and lower boundaries are 
both solid surfaces, the two values of F wake at a particu- 
lar streamwise station are combined using the following 
averaging formula: 


( F w*ke^ + ^ F wakc^2^ 2 


wake 


f,+f. 


(A. 10) 


The quantities (F wlke ), and (F wtke ) 2 are the separate 
values computed at the lower and upper boundaries 
using equation (A.6). The functions fj and {2 are defined 
by 


wake 


v F 

Jmax 1 


max max 


p v* 1- 

c wk v diffp 


for wall bounded flows 

(A.6) 

for free turbulent flows 


where the constant C wk is 0.25 and 

V diff = |vL„-|vLi„ (A.7) 

where V is the total velocity vector. The quantity F max 
is the maximum value of 


F(y) = 


yM 




1 -e 


for wall bounded flows 

(A.8) 


yM 


for free turbulent flows 


and y m)UC is the value of y corresponding to F max- F Kleb 
is the Klebanoff intermittency factor which accounts for 
the experimentally observed phenomenon that as the free 
stream is approached, the fraction of time the flow is 
turbulent decreases. It is given by 



where is a constant taken as 0.3. 



(A. 11) 


The constant n is set equal to 2.0, y\ and y 2 are the 
normal distances to the bottom and top surfaces, respec- 
tively, and Dj and D 2 are the normal distances from the 
two surfaces to the location of |v| max . In addition, the y/ 
y max value used in equation (A.9) for F^^ is computed 
for both surfaces and the minimum value is used. These 
values ofF K i cb and F wakc are then used in equation (A.5) 
to compute (H t ) outcf . 

A ppendix B: Artificial Viscosity and Grid Convergence 

High frequency nonlinear instabilities can appear 
as the Proteus solution develops. For example, physical 
phenomena, such as shock waves, can cause instabilities 
when they are captured by the finite difference 
algorithm. In addition, high Reynolds number flows 
may have oscillations resulting from the odd-even 
decoupling inherent in the use of central spatial differ- 
encing of the convective terms. Artificial viscosity may 
be used to suppress these oscillations. The two artificial 
viscosity models in Proteus are the constant coefficient 
model of Steger 29 and the nonlinear coefficient model of 
Jameson et al. 28 The implementation of these models in 
generalized nonorthogonal coordinates was taken from 
Pulliam. 31 


3. Multiple Boundaries 

If both walls in a given coordinate direction are 
solid surfaces, the turbulence model equations are 
applied separately at each surface and then averaged. 
The two outer regions overlap, and it assumed that the 
two inner regions do not overlap. The averaging proce- 
dure deals with the F wakc function. For example, in the 


1 . Constant Coefficient Model 

The constant coefficient model uses a combina- 
tion of explicit and implicit smoothing. The standard 
explicit artificial viscosity uses fourth-order differ- 
ences. Second-order explicit artificial viscosity, which 
provides more smoothing, is also available in Proteus, 
however it is rarely used because it introduces a large 
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error. The implicit smoothing is second order and is 
used to extend the linear stability bound of the fourth- 
order explicit smoothing. 

The explicit artificial viscosity is implemented in 
the Proteus alternating direction implicit (ADI) 
algorithm 2 by adding the following terms to the right- 
hand side source term for the first ADI sweep. (The 
governing equations of Proteus are given in detail in 
reference 1 .) 

e< 2) * At 

-V-(W + w>- 

f4) <B.l) 

4 4) a t 

-^j-UV^Q+^A^Q] 

where e^ 2) and e^ 4) are the second- and fourth-order 
explicit artificial viscosity coefficients, £ and r\ are the 
computational coordinate directions, Q is the vector of 
dependent variables and J is the Jacobian of the coordi- 
nate transformation. The symbols V and A are the 
standard backward and forward first difference opera- 
tors. 

The implicit artificial viscosity is implemented 
by adding the following terms to the left-hand side of the 
governing equation. 


♦v5.,*5Ks 

The expression y is given by 


^Q-e^A^V^Q).} 
Q - e* 4) A V A Q) } 


(B3) 


V = V* + V y - (B4) 

where v and y are spectral radii defined by 

x y 


|U|+a^ + ^ 

= Al 

ivi+a7n 2 +n 2 

Arj 


(B.5) 


The second- and fourth- order nonlinear artificial 
viscosity coefficients are a function of the pressure field. 
In the £ direction, they are given by 


(e^ 2) ) t = KjAt maxfa^j.Cj.o.,,) (B.6a) 

(e| 4) ) t = max[0,K 4 At- (e^ 2) ) J (B.6b) 


zA t ^ # _ 

- j [V,<"G>] 

(B.2a) 


(B.2b) 


Equation (B.2a) is added for the first ADI sweep and 
equation (B.2b) is added for the second ADI sweep. The 
constant is the implicit artificial viscosity coefficient. 

The optimum values of the coefficients e£ 2) , e^ 4) 
and vary from problem to problem. They should be 
small so as not to corrupt the physical solution, yet large 
enough to damp any instabilities. The Proteus User’s 
Guide 1 recommends starting values of e*. 4) =1.0, 
ef =1.0 and e 7 =2.0. 


where 


G i^ 


Pm~ 2 Pi + Pi-i 

Pi + t + 2 Pi + Pi_, 


(B.7) 


Similar formulas are used in the q direction. 

The parameter a is a pressure gradient scaling 
parameter which increases the amount of second-order 
smoothing relative to fourth-order near shock waves. 
The parameters k 2 and 1 C 4 are user-specified constants. 
As with the constant coefficient model, the optimum 
values of k 2 and K 4 are problem dependent Typical 
values range from K4=0.005 and k 2 =0.01 for flows with 
no shocks, to K4=0.0004 and k 2 =0.1 for flows with 
shocks . 1 Pulliam gives k 2 =0.25 and K 4 = 0.01 as typical 
values for an Euler analysis . 1 * 31 


2. Nonlinear Coefficient Model 

The nonlinear coefficient artificial viscosity is 

explicit and contains second and fourth-order differ- 
ences. The following terms are added to the right-hand 
side of the governing equations. 


3. Comments on Artificial Viscosity 

As previously mentioned, artificial viscosity is 
generally used to minimize oscillations which occur 
when computing high Reynolds number flows and flows 
with shock waves. The artificial viscosity coefficients 
should be as small as possible so as not to corrupt the 
solution, yet large enough to damp the nonphysical insta- 
bilities. Optimum values of the artificial viscosity 
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coefficients vary from problem to problem; the coeffi- 
cients used to calculate the flows presented herein were 
selected based on values used for similar cases, as given 
in the Proteus User’s Manual. 1 Some representative test 
cases were evaluated to insure that the chosen artificial 
viscosity did not corrupt the physical characteristics of 
the flow. 

The flat plate flows were run using the constant 
coefficient model with ££^=1.0, ££^=0.0 and £/=2,0. 
For these flows, it was possible to run Proteus with zero 
artificial viscosity, however the solutions took two to 
four times longer to converge. Upon close examination, 
these solutions did not agree as closely with experimen- 
tal data as the solutions computed using artificial viscos- 

ity- 

For the diffuser flows, the artificial viscosity 
effects were examined for the weak shock case. As 
mentioned earlier, the nonlinear coefficient model was 
used. A value of *2=0.1 was used for the entire calcula- 
tion, with *4=0.005 while the exit pressure was changing 
(i.e, for the first 6000 iterations), and k 4 =0.0004 for the 
remaining 3000 iterations, which were at a constant exit 
pressure. It was not possible to compute this flow 
without artificial viscosity, so the effect of doubling and 
halving the smoothing coefficients was examined. The 
static pressure distribution on the top and bottom walls 
for this comparison (computed using MMLPG) is given 
in figure 24. The solution computed using half of the 
original artificial viscosity was nearly identical to the 
original solution, indicating that the originally chosen 
artificial viscosity is reasonable for this flow. Doubling 
the artificial viscosity gave a less desirable result in that 
the normal shock was smeared over a greater number of 
grid points. 

4. Grid Convergence 

Grid convergence is an important factor in the 
accuracy of a CFD calculation. The grids used to make 

the flat plate and transonic diffuser calculations were 
assessed to insure their grid independence. For the zero 
pressure gradient flat plate calculations, a 101x101 grid 
was initially chosen. The size of this grid was systemat- 
ically reduced in each direction in order to find the 
coarsest grid that would give a solution which would not 
change if additional grid points were added. The 51x51 
grid shown in figure 5 was chosen based on this proce- 
dure. 

The grid used to make the transonic diffuser 
calculations had been used previously by others,^ 2 ^ so it 
is probable that this grid gives a reliable solution. As an 
added check, the number of grid points in each direction 
was doubled, and the resulting 162x101 grid was used to 


compute the no shock flow using MMLPG. A compari- 
son of these results with the results obtained using the 
81x51 grid of figure 17 is given in figure 25 and 
indicates that the 81x51 grid is reliable. 


Table 1. Parameters used in pressure gradient 
modifications. 


Pressure 

Gradient 

Strength 

p 

g 3 

g 5 

°6 

zero 

0 

1.00 

23,300 

-7.75xl0 6 

mild 

1 

1.25 

30,100 

-1.16xl0 7 

strong 

5 

1.53 

33,800 

-2.09x1 0 7 


Table 2. Computational times for flat plate flows. 


(a) Zero Pressure Gradient 


Model 

Iterations 

sec./iter./grid 

point 

BLM 

2000 

2.02x1 0^ 

MMLPG 

2000 

2.00x1 O' 4 5 


(b) Mild Pressure Gradient 


Model 

Iterations 

scc./iter./grid 

point 

BLM 

3000 

2.14xl0" 5 

MMLPG 

3000 

2.17zl0' 5 


(c) Strong Pressure Gradient 


Model 

Iterations 

sec./iter./grid 

point 

BLM 

18,000 

2.14xl0‘ 5 

MMLPG 

10,000 

1.96xl0‘ 5 


Table 3. Shock location and Mach number, weak 
shock case. 


Turbulence 

Model 

Shock Mach 
Number 

Shock Loca- 
tion (x/H) 

MMLPG 

1.233 

1.57 

BLM 

1.228 

1.49 

Experiment 22 

1.235 

1.41 
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Table 4. Maximum Mach number, no 
shock case 


Turbulence 

Model 

Maximum 
Mach Number 

MMLPG 

0.881 

BLM 

0.976 


Table 5. Shock location and Mach number, strong 
shock case. 


Turbulence 

Model 

Shock Mach 
Number 

Shock Loca- 
tion (x/H) 

MMLPG 

1.626 

3.13 

BLM 

1.665 

2.90 

Experiment 22 

1.353 

1.98 









y/8 


Figure 6. Velocity-defect profiles for zero pressure gradient 
flow, Re„=7x10 . 



(a) MML (b) BLM 


Figure 7. Velocity-defect profiles for zero pressure flow at 
three Reynolds numbers. 



(a) MML (b) BLM 

Figure 8. Turbulent viscosity for zero pressure gradient flat 
plate flow at three Reynolds numbers. 
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0.0050 



Figure 9. Velocity-defect profiles for zero pressure gradient 
flow, MML of equations (9) and (11). 
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Figure 10. Local skin friction coefficient for zero pressure 
gradient flow: MML of equations (9) and (11 ). 




y/5 

Figure 12. Velocity-defect for zero pressure gradient flow; 
MML of equations (13) through (17). 



(a) MML of equations (13) through (17) (b) BLM 

Figure 13. Velocity-defect for mild pressure gradient flow. 
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(a) MML of equations (13) through (17) 

Figure 14. Velocity-defect for strong pressure gradient flow. 
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(b) BLM 
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(a) Zero pressure gradient 
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(b) Mild pressure gradient 
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(c) Strong pressure gradient 

Figure 15. Velocity-defect profiles computed using MMLPG. 
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(a) Top wall (b) Bottom wall 

Figure 24. Static pressure, computed using MMLPG and 
three different amounts of artificial viscosity. 


20 
















5 . 


(a) Top wall (b) Bottom wall 

Figure 25. Static pressure computed using MMLPG and two 
different grids; no shock case. 












REPORT DOCUMENTATION PAGE 


Form Approved 
OMB No. 0704-0188 


■ p - ubUc ^porting burden tor this collection of info rmation is estimated to average 1 hour per response, including the time for wjj^na Instny Ions searching existing dal^ources, 
aatherino wwj^maintalning the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other 
collection of information including suggestions for reducing this burden, to Washington Headquarters Servioes, Directorate tor Infor^tion C^erations and Reports, 1 2 1 5Jefferson 
Davis Highway, Suite 1204, Arlington? V A 22202-4302. and to the Office of Management and Budget. Paperwork Reduction Project (0704-0188), Washington, DC 20503. 

i. AGENCY USE ONLY (Leave blank) 1 2. REPORT DATE 1 3. REPORT TYPE AND DATES COVERED ~~ ~ 


November 1994 


4. TITLE AND SUBTITLE 

A Modified Mixing Length Turbulence Model for Zero and 
Adverse Pressure Gradients 


6. AUTHOR(S) 

J.M. Conley and B.P. Leonard 


Technical Memorandum 


5. FUNDING NUMBERS 


WU- 505-62-52 


7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 

National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 44135-3191 


8. PERFORMING ORGANIZATION 
REPORT NUMBER 


E-9220 


9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESSEES) 

National Aeronautics and Space Administration 
Washington, D.C. 20546-0001 


10. SPONSORING/MONITORING 
AGENCY REPORT NUMBER 


NASA TM- 106772 


Prepared for the 30th Joint Propulsion Conference cosponsored by AIAA, ASME, SAE, and ASEE, Indianapolis, Indiana, June 27-29, 1994 under the 
title “Modification of the MML Turbulence Model for Adverse Pressure Gradient Flows.” J.M. Conley, NASA Lewis Research Center, and B.P. Leonard, 
The University of Akron, Akron, Ohio 44325. This material is a portion of NAS ATM- 106544 which was submitted by J.M. Conley as her Masters 
t hesis. Responsible person, J.M. Conley, organization code 2660, (216) 433—2188. 

12a. DISTRIBUTION/ AVAILABILITY STATEMENT ‘ ” 12b. DISTRIBUTION CODE 

Unclassified -Unlimited 
Subject Category 34 

13. ABSTRACT (Maximum 200 words) 

The modified mixing length (MML) turbulence model was installed in the Proteus Navier-Stokes code, then modified to 
make it applicable to a wider range of flows typical of aerospace propulsion applications. The modifications are based 
on experimental data for three flat-plate flows having zero, mild adverse, and strong adverse pressure gradients. Three 
transonic diffuser test cases were run with the new version of the model in order to evaluate its performance. All results 
are compared with experimental data and show improvements over calculations made using the Baldwin-Lomax 
turbulence model, the standard algebraic model in Proteus. 


14. SUBJECT TERMS 

Navier-Stokes; Turbulent boundary layer; Two-dimensional flow; Mixing length 
flow theory; Pressure gradients; Flat plates; Computational fluid dynamics 


17. SECURITY CLASSIFICATION 
OF REPORT 

Unclassified 

NSN 7540-01-280-5500 


18. SECURITY CLASSIFICATION 
OF THIS PAGE 

Unclassified 


19. SECURITY CLASSIFICATION 
OF ABSTRACT 

Unclassified 


15. NUMBER OF PAGES 

23 


16. PRICE CODE 

A03 


, 20. LIMITATION OF ABSTRACT 


Standard Form 298 (Rev. 2-89) 
Proscribed by ANSI Std. Z39-18 
298-102 







